<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/1999/REC-html401-19991224/loose.dtd">
<html lang="en">
<head>
	<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
	<title>Simulations In Three Dimensions Example (k-Wave)</title>
	<link rel="stylesheet" href="kwavehelpstyle.css" type="text/css">
	<meta name="description" content="Simulations In Three Dimensions Example.">
</head>

<body><div class="content">

<h1>Simulations In Three Dimensions Example</h1>

<p>This example provides a simple demonstration of using k-Wave for the simulation and detection of the pressure field generated by an initial pressure distribution within a three-dimensional heterogeneous propagation medium. It builds on the <a href="example_ivp_homogeneous_medium.html">Homogeneous Propagation Medium</a> and <a href="example_ivp_heterogeneous_medium.html">Heterogeneous Propagation Medium</a> examples.</p>

<p>
    <ul>
        <li><a href="matlab:edit([getkWavePath('examples') 'example_ivp_3D_simulation.m']);" target="_top">Open the file in the MATLAB Editor</a></li>
        <li><a href="matlab:run([getkWavePath('examples') 'example_ivp_3D_simulation']);" target="_top">Run the file in MATLAB</a></li>
    </ul>
</p>

<h2>Contents</h2>
<div>
	<ul>
        <li><a href="#heading2">Creating the k-space grid and defining the medium properties</a></li>
        <li><a href="#heading3">Defining the initial pressure distribution</a></li>
        <li><a href="#heading4">Defining the sensor mask</a></li>
        <li><a href="#heading5">Running the simulation</a></li>
        <li><a href="#heading6">A note on Cartesian sensor masks</a></li>
    </ul>
</div>

<a name="heading2"></a>
<h2>Creating the k-space grid and defining the medium properties</h2>

<p>Simulations in three-dimensions are performed in an analogous fashion to those in one and two dimensions. The medium discretisation is again performed by <code><a href="kWaveGrid.html">kWaveGrid</a></code> with two additional inputs for the z-dimension. Similarly, matrices for the medium properties also have an extra dimension. The matrix indexing convention (x, y, z) is an extension to that used in two-dimensions (x, y) where matrix row and column data correspond to the x and y axes, respectively.</p>

<pre class="codeinput">
<span class="comment">% create the computational grid</span>
Nx = 64;            <span class="comment">% number of grid points in the x direction</span>
Ny = 64;            <span class="comment">% number of grid points in the y direction</span>
Nz = 64;            <span class="comment">% number of grid points in the z direction</span>
dx = 0.1e-3;        <span class="comment">% grid point spacing in the x direction [m]</span>
dy = 0.1e-3;        <span class="comment">% grid point spacing in the y direction [m]</span>
dz = 0.1e-3;        <span class="comment">% grid point spacing in the z direction [m]</span>
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);

<span class="comment">% define the properties of the propagation medium</span>
medium.sound_speed = 1500 * ones(Nx, Ny, Nz);	<span class="comment">% [m/s]</span>
medium.sound_speed(1:Nx/2, :, :) = 1800;        <span class="comment">% [m/s]</span>
medium.density = 1000 * ones(Nx, Ny, Nz);       <span class="comment">% [kg/m^3]</span>
medium.density(:, Ny/4:end, :) = 1200;          <span class="comment">% [kg/m^3]</span>
</pre>

<a name="heading3"></a>
<h2>Defining the initial pressure distribution</h2>

<p>As in one and two dimensions, the initial pressure distribution is defined as a matrix of arbitrary numbers the same size as the computational grid. Here <code><a href="makeBall.html">makeBall</a></code> is used to create an initial pressure distribution of two small filled balls with different diameters.</p>

<pre class="codeinput">
<span class="comment">% create initial pressure distribution using makeBall</span>
ball_magnitude = 10;    <span class="comment">% [Pa]</span>
ball_x_pos = 38;        <span class="comment">% [grid points]</span>
ball_y_pos = 32;        <span class="comment">% [grid points]</span>
ball_z_pos = 32;        <span class="comment">% [grid points]</span>
ball_radius = 5;        <span class="comment">% [grid points]</span>
ball_1 = ball_magnitude * makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

ball_magnitude = 10;    <span class="comment">% [Pa]</span>
ball_x_pos = 20;        <span class="comment">% [grid points]</span>
ball_y_pos = 20;        <span class="comment">% [grid points]</span>
ball_z_pos = 20;        <span class="comment">% [grid points]</span>
ball_radius = 3;        <span class="comment">% [grid points]</span>
ball_2 = ball_magnitude * makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

source.p0 = ball_1 + ball_2;
</pre>

<a name="heading4"></a>
<h2>Defining the sensor mask</h2>

<p>Again the sensor mask can be given as a list of Cartesian coordinates, a binary mask, or the grid coordinates of two opposing corners of a cuboid. Several functions for creating three-dimensional shapes are included in the k-Wave toolbox. For example, <code><a href="makeSphere.html">makeSphere</a></code> which returns a binary map of single pixel spherical shell created using the midpoint circle algorithm, and <code><a href="makeCartSphere.html">makeCartSphere</a></code> which returns the Cartesian location of an evenly distributed array of an arbitrary number of points on the surface of a sphere. In this example, a linear sensor array is explicitly created.</p>

<pre class="codeinput">
<span class="comment">% define a series of Cartesian points to collect the data</span>
x = (-22:2:22) * dx;          <span class="comment">% [m]</span>
y = 22 * dy * ones(size(x));    <span class="comment">% [m]</span>
z = (-22:2:22) * dz;          <span class="comment">% [m]</span>
sensor.mask = [x; y; z];
</pre>

<p>A visualisation of the initial pressure distribution and the sensor mask using <code><a href="cart2grid.html">cart2grid</a></code> and <code><a href="voxelPlot.html">voxelPlot</a></code> is given below.</p>

<img vspace="5" hspace="5" src="images/example_ivp_3D_simulation_01.png" style="width:560px;height:420px;" alt="">

<a name="heading5"></a>
<h2>Running the simulation</h2>

<p>The computation is invoked by calling <code><a href="kspaceFirstOrder3D.html">kspaceFirstOrder3D</a></code> with the inputs defined above. By default, a visualisation of the wave-field in the horizontal, median, and frontal planes are displayed, along with a status bar.</p>

<pre class="codeinput">
<span class="comment">% input arguments</span>
input_args = {'PlotLayout', true, 'PlotPML', false, ...
    'DataCast', 'single', 'CartInterp', 'nearest'};

<span class="comment">% run the simulation</span>
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
</pre>

<p>By setting the optional input parameter <code>'PlotLayout'</code> to <code>true</code>, a plot of the initial pressure, and sound speed and density (if heterogenous) are produced. To remove the PML from the display, the optional input <code>'PlotPML'</code> is also set to <code>false</code>.</p>

<img vspace="5" hspace="5" src="images/example_ivp_3D_simulation_02.png" style="width:560px;height:420px;" alt="">
<img vspace="5" hspace="5" src="images/example_ivp_3D_simulation_03.png" style="width:560px;height:420px;" alt="">
<img vspace="5" hspace="5" src="images/example_ivp_3D_simulation_04.png" style="width:560px;height:420px;" alt="">

<p>As the function runs, status updates and computational parameters are printed to the command line. By setting the optional input <code>'DataCast'</code> to <code>'single'</code>, the computational time is decreased (see the <a href="example_na_optimising_performance.html">Optimising k-Wave Performance Example</a> for more details).</p>

<pre class="codeinput">
Running k-Wave simulation...
  start time: 13-Apr-2017 22:02:40
  reference sound speed: 1800m/s
  dt: 16.6667ns, t_end: 7.3833us, time steps: 444
  input grid size: 64 by 64 by 64 grid points (6.4 by 6.4 by 6.4mm)
  maximum supported frequency: 7.5MHz
  smoothing p0 distribution...
  casting variables to single type...
  precomputation completed in 2.7428s
  starting time loop...
  estimated simulation time 15.6821s...
  simulation completed in 16.4262s
  reordering Cartesian measurement data...
  total computation time 19.286s
</pre>

<p>When the time loop has completed, the function returns the recorded time series at each of sensor locations defined by <code>sensor.mask</code>. The ordering is again dependent on whether a Cartesian or binary sensor mask is used. A visualisation is given below.</p>

<img vspace="5" hspace="5" src="images/example_ivp_3D_simulation_05.png" style="width:560px;height:420px;" alt="">

<a name="heading6"></a>
<h2>A note on Cartesian sensor masks</h2>

<p>If a Cartesian sensor mask is used, the values of the acoustic field at the sensor points are obtained at each time step using interpolation. The interpolation mode can be changed by setting the optional input parameter <code>'CartInterp'</code> to <code>'linear'</code> (the default) or <code>'nearest'</code> (used in this example). During the simulation, there is only a small performance difference between using Cartesian and binary sensor masks. However, if a Cartesian sensor mask with linear interpolation is used, the calculation of the triangulation points can significantly lengthen the precomputation time, particularly for 3D simulations (see the output below for the same example run with <code>'CartInterp'</code> set to <code>'linear'</code>). A Cartesian sensor mask can be converted to a binary sensor mask (and vice versa) using the functions <code><a href="cart2grid.html">cart2grid</a></code> and <code><a href="grid2cart.html">grid2cart</a></code>.</p>

<pre class="codeinput">
Running k-Wave simulation...
  start time: 18-Oct-2011 10:08:44
  reference sound speed: 1800m/s
  dt: 16.6667ns, t_end: 6.15us, time steps: 370
  input grid size: 64 by 64 by 64 grid points (6.4 by 6.4 by 6.4mm)
  maximum supported frequency: 7.5MHz
  smoothing p0 distribution...
  calculating Delaunay triangulation...
  precomputation completed in 25.1272s
  starting time loop...
  estimated simulation time 44.8s...
  simulation completed in 50.4407s
  total computation time 1min 15.598s
</pre>

<p>In both two and three dimensional simulations, if the optional input <code>'CartInterp'</code> is set to <code>'nearest'</code>, the Cartesian points are internally mapped onto a binary grid. If more than one Cartesian sensor point maps to the same location, the duplicated points are discarded. Consequently, if the Cartesian sensor mask has a dense array of points, the number of returned time-series may not correspond to the original number of Cartesian sensor points. Details of discarded duplicate points are printed to the command line. This situation can be avoided by increasing the resolution of the computational grid, or by reducing the packing of the Cartesian sensor points.</p>

</div></body></html>